/**********************************************************************
int select_matrix (int N, int matrix[N][N], int order[N])
This function is wrote for solving the following optimization problem:
Give a N x N matrix and limit only one element could be selected in each
row and each column, how to select N elements to let the sum of them reach
the minimum? The return value gives the minimal sum and the selection is
matrix[order[i]][i] (i form 0 to N - 1).
**********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

struct elimate {
	int diff;
	int idx;
};

/* simplify the select matrix : if in i-th column only the j-th element is selected, then the other elements in the j-th is elimated */
int elimate2 (int N, int select[N + 1][N + 1], int order[N])
{
	int i, j, pos, counter, flag;

	counter = 0;

	flag = 1;
	while (flag)
	{
		flag = 0;

		for (i = 0; i < N; i++)
		{
			if (select[i][N] == 1)
			{
				for (j = 0; j < N; j++)
				{
					if (select[i][j] == 1)
					{
						order[j] = i;
						pos = j;
						break;
					}
				}
        
				if (select[N][pos] > 1)
				{
					counter += select[N][pos] - 1;
					select[N][N] -= select[N][pos] - 1;
					select[N][pos] = 1;
        
					for (j = 0; j < N; j++)
					{
						if (j == i)
							continue;
        
						if (select[j][pos] == 1)
						{
							select[j][pos] = 0;
							select[j][N] --;
							flag = 1;
						}
					}
				}
			}
		}
        
		for (i = 0; i < N; i++)
		{
			if (select[N][i] == 1)
			{
				for (j = 0; j < N; j++)
				{
					if (select[j][i] == 1)
					{
						order[i] = j;
						pos = j;
						break;
					}
				}
        
				if (select[pos][N] > 1)
				{
					counter += select[pos][N] - 1;
					select[N][N] -= select[pos][N] - 1;
					select[pos][N] = 1;
        
					for (j = 0; j < N; j++)
					{
						if (j == i)
							continue;
        
						if (select[pos][j] == 1)
						{
							select[pos][j] = 0;
							select[N][j] --;
							flag = 1;
						}
					}
				}
			}
		}
	}

	return counter;
}

int elimate1 (int N, struct elimate eli_col[N][N][N], struct elimate eli_row[N][N][N], int select[N + 1][N + 1], int order[N])
{
	int i, j, k, counter, flag;

	counter = 0;

	for (i = 0; i < N; i++)
	{
		for (j = 0; j < N; j++)
		{
			if (i == j)
				continue;

			if (select[N][i] >= 2)
			{
				flag = 0;
				for (k = 0; k < N - 1; k++)
				{
					if (select[eli_col[i][j][k].idx][j] == 0)
					{
						continue;
					}
					else
					{
						if (select[eli_col[i][j][k].idx][i] == 1)
						{
							if (eli_col[i][j][k].diff > eli_col[i][j][k + 1].diff)
								flag = 1;
						}
						break;
					}
				}

				/* elimate matrix[eli_col[i][j][k].idx][i] (because matrix[eli_col[i][j][k].idx][j] is always a better choice) */
				if (flag == 1)
				{
					select[eli_col[i][j][k].idx][i] = 0;
					select[N][i] --;
					select[eli_col[i][j][k].idx][N] --;
					select[N][N] --;
					counter ++;
					if (select[N][i] == 1 || select[eli_col[i][j][k].idx][N] == 1)
						counter += elimate2 (N, select, order);
				}
			}

			if (select[i][N] >= 2)
			{
				flag = 0;
				for (k = 0; k < N - 1; k++)
				{
					if (select[j][eli_row[i][j][k].idx] == 0)
					{
						continue;
					}
					else
					{
						if (select[i][eli_row[i][j][k].idx] == 1)
						{
							if (eli_row[i][j][k].diff > eli_row[i][j][k + 1].diff)
								flag = 1;
						}
						break;
					}
				}

				/* elimate matrix[i][eli_row[i][j][k].idx] (because matrix[j][eli_row[i][j][k].idx] is always a better choice) */
				if (flag == 1)
				{
					select[i][eli_row[i][j][k].idx] = 0;
					select[N][eli_row[i][j][k].idx] --;
					select[i][N] --;
					select[N][N] --;
					counter ++;
					if (select[i][N] == 1 || select[N][eli_row[i][j][k].idx] == 1)
						counter += elimate2 (N, select, order);
				}
			}
		}
	}

	return counter;
}

int compare (struct elimate *a, struct elimate *b)
{
	if (a->diff > b->diff)
		return -1;
	else
		return 1;
}

void setup_elimate (int N, int matrix[N][N], struct elimate eli_col[N][N][N], struct elimate eli_row[N][N][N])
{
	int i, j, k;

	for (i = 0; i < N; i++)
	{
		for (j = 0; j < N; j++)
		{
			if (i == j)
				continue;

			for (k = 0; k < N; k++)
			{
				eli_col[i][j][k].diff = matrix[k][i] - matrix[k][j];
				eli_col[i][j][k].idx = k;

				eli_row[i][j][k].diff = matrix[i][k] - matrix[j][k];
				eli_row[i][j][k].idx = k;
			}

			qsort (eli_col[i][j], N, sizeof (struct elimate), compare);
			qsort (eli_row[i][j], N, sizeof (struct elimate), compare);
		}
	}
}

int calc_score (int N, int matrix[N][N], int order[N])
{
	int i;
	int score;

	score = 0.0;
	for (i = 0; i < N; i++)
	{
		score += matrix[order[i]][i];
	}

	return score;
}

int elimate (int N, int matrix[N][N], struct elimate eli_col[N][N][N], struct elimate eli_row[N][N][N], int select[N + 1][N + 1], int order[N])
{
	int i, j, k, counter, flag;
	int score;

	score = 0.0;

	for (;;)
	{
		counter = elimate1 (N, eli_col, eli_row, select, order);
		if (counter == 0 || select[N][N] == N)
			break;
	}

	flag = 0;
	for (i = 0; i < N; i++)
	{
		if (select[i][N] * select[N][i] == 0)
		{
			flag = 0;
			break;
		}
		else
		{
			flag += select[i][N] * select[N][i];
		}
	}

	if (flag == N)
	{
		score = calc_score (N, matrix, order);
	}
	else if (flag == 0)
	{
		score = -1.0;
	}
	else
	{
		int select2[N + 1][N + 1], order2[N];
		int select_min[N + 1][N + 1], order_min[N];
		int score2, score_min;

		for (i = 0; i < N; i++)
		{
			if (order[i] == -1)
				break;
		}

		counter = 0;
		for (j = 0; j < N; j++)
		{
			if (select[j][i] == 0)
				continue;

			memcpy (select2, select, sizeof (int) * (N + 1) * (N + 1));
			memcpy (order2, order, sizeof (int) * N);

			order2[i] = j;
			select2[N][N] -= select2[N][i] - 1;
			select2[N][i] = 1;
			for (k = 0; k < N; k++)
			{
				if (k == j)
					continue;

				if (select2[k][i] == 1)
				{
					select2[k][i] = 0;
					select2[k][N] --;
				}
			}
			elimate2 (N, select2, order2);

			flag = 0;
			for (k = 0; k < N; k++)
			{
				if (select[k][N] * select[N][k] == 0)
				{
					flag = 0;
					break;
				}
				else
				{
					flag += select[k][N] * select[N][k];
				}
			}

			if (flag == 0)
			{
				continue;
			}
			else if (flag == N)
			{
				score2 = calc_score (N, matrix, order2);
			}
			else
			{
				score2 = elimate (N, matrix, eli_col, eli_row, select2, order2);
			}

			if (score2 < 0)
				continue;

			counter ++;

			if (counter == 1)
			{
				memcpy (select_min, select2, sizeof (int) * (N + 1) * (N + 1));
				memcpy (order_min, order2, sizeof (int) * N);
				score_min = score2;
			}
			else if (score2 < score_min)
			{
				memcpy (select_min, select2, sizeof (int) * (N + 1) * (N + 1));
				memcpy (order_min, order2, sizeof (int) * N);
				score_min = score2;
			}
		}

		if (counter > 0)
		{
			memcpy (select, select_min, sizeof (int) * (N + 1) * (N + 1));
			memcpy (order, order_min, sizeof (int) * N);
			score = score_min;
		}
		else
		{
			score = -1.0;
		}
	}

	return score;
}

int select_matrix (int N, int matrix[N][N], int order[N])
{
	int score;
	int i, j, counter;
	/* select[i][j] : 1 means that matrix[i][j] is a candidate; 0 means that matrix[i][j] is elimated; select[N][i]: sum of select[j][i] (j: from 0 to N-1); select[i][N]: sum of select[i][j] (j: from 0 to N-1); select[N][N]: sum of select[i][j] (i,j: from 0 to N-1) */
	int select[N + 1][N + 1];
	/* eli_col[i][j][k] : the difference when swapping i-th column and j-th column, sorted */
	struct elimate eli_col[N][N][N];
	/* eli_row[i][j][k] : the difference when swapping i-th row and j-th row, sorted */
	struct elimate eli_row[N][N][N];

	setup_elimate (N, matrix, eli_col, eli_row);

	for (i = 0; i < N; i++)
	{
		order[i] = -1;
	}

	for (i = 0; i < N; i++)
	{
		for (j = 0; j < N; j++)
		{
			select[i][j] = 1;
		}
	}
	for (i = 0; i < N + 1; i++)
	{
		select[N][i] = N;
		select[i][N] = N;
	}
	select[N][N] = N * N;

	score = elimate (N, matrix, eli_col, eli_row, select, order);

	return score;
}
